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1 .0  IDENTIFICATION  OF  THE  PROBLEM 

The  potential  health  effects  of  exposure  to  sources  of  non-ionizing  radio  frequency  (RF)  electro¬ 
magnetic  energy  are  an  area  of  continuing  interest  within  the  military  and  civilian  communities. 
Common  sources  of  RF  exposure  include  mobile  and  cellular  telephones,  magnetic  resonance 
imaging  (MRI)  systems,  wireless  local  area  networks,  transmitting  antennae,  and  other  civilian 
and  military  communications  and  radar  systems.  One  established  biological  effect  of  RF 
exposure  is  tissue  heating;  temperature  increases  as  little  as  4°C  above  normal  body  temperature 
can  have  potentially  devastating  effects  on  living  tissue  [31].  This  report  discusses  the 
development  of  a  prototype  of  an  anatomically  realistic  thermal  computer  code  capable  of 
predicting  the  thermal  effects  of  known  RF  sources  on  human  and  animal  tissue. 

The  Institute  of  Electrical  and  Electronics  Engineers  (IEEE)  recently  adopted  revised  standards 
for  human  exposure  to  RF  fields,  specifying  maximum  permissible  exposure  (MPE)  levels.  The 
MPE  level,  for  a  given  exposure  frequency,  is  stated  in  terms  of  the  specific  absorption  rate 
(SAR)  -  the  rate  of  energy  absorption  per  unit  mass  of  tissue.  The  MPE  depends  upon  the 
duration  of  exposure  and  whether  the  exposure  is  “controlled”  to  avoid  physiological  hazards 
[28].  The  rate  of  temperature  increase  in  tissue  is  a  function  of  the  specific  absorption  rate.  The 
local  rate  of  temperature  increase  in  a  homogeneous  phantom  model  is  approximately  linear  with 
the  SAR,  for  short-duration  exposures  [27].  Phantom  models  however,  are  of  little  value  in 
predicting  tissue  heating  effects  of  RF  exposure  in  humans  or  animals.  The  effects  of  metabolic 
heating,  blood  flow,  and  tissue  heterogeneity  complicate  the  relationship  between  the  rate  of 
temperature  increase  and  the  local  SAR,  such  that  a  simple  linear  model  is  not  generally 
appropriate. 

The  anatomical  complexity  and  high  degree  of  spatial  resolution  required  to  capture  predicted 
SAR  “hot  spots”  have  led  to  the  development  of  realistic  anatomical  data  sets  [1,2].  These  voxel- 
based  descriptions  are  created  by  dividing  the  space  occupied  by  the  object  being  described  (e.g., 
the  human  body)  into  a  three  dimensional  grid  of  small,  equal-sized  volume  elements  known  as 
voxels.  By  definition,  each  voxel  can  contain  only  one  type  of  tissue;  a  very  large  number  of 
small  voxels  (e.g.,  cubes  with  1mm  edges)  is  therefore  needed  to  capture  enough  detail  to 
accurately  predict  SARs.  The  properties  of  the  tissue  in  each  voxel  and  the  predicted  SARs 
constitute  the  input  to  a  thermal  model,  which  then  predicts  tissue  temperatures.  Consequently, 
general-purpose  off-the-shelf  thermal  codes  are  quickly  overwhelmed  by  the  size  of  thermal 
problems  that  are  derived  from  voxel-based  anatomical  descriptions.  ThermoAnalytics  proposed 
to  develop  new  thermal  software  to  exploit  the  simple  structure  of  the  voxel-based  description. 
The  resulting  thermal  code  is  compact,  accurate,  and  able  to  accommodate  extremely  large  data 
sets  as  input. 

This  report  addresses  the  issue  of  thermoregulation  and  bio-heat  transfer  modeling  as  it  pertains 
to  complex  voxel-based  anatomical  descriptions  with  thermal  loads  from  electromagnetic  (EM) 
irradiation.  The  level  and  location  of  heat  loading  are  derived  using  the  finite  difference  time 
domain  (FDTD)  method  and  applied  to  the  model  as  a  specific  absorption  rate.  The  applicability 
of  the  resulting  computer  code  is  not  limited  to  RF  heating  however,  but  extends  to  a  variety  of 
environmental  and  therapeutic  heating  and  cooling  applications. 
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1.1  Background 

The  development  of  accurate,  whole-body,  predictive  models  of  heat  transfer  in  humans  and 
primates  would  be  an  extremely  important  accomplishment,  both  scientifically  and  for  the 
potential  economic  benefit  deriving  from  the  use  of  the  models.  While  whole-body  models  have 
been  described  in  the  literature  [2,10],  the  lack  of  anatomical  detail  and  precision,  and  limited 
thermoregulatory  mechanisms,  means  their  applications  are  limited  to  describing  relatively  gross 
thermal  response. 


Thermoregulation  is  an  ensemble  of  physiological  processes  that  differ  among  species,  and 
which  collectively  complicate  the  development  of  predictive  models.  Conduction,  convection, 
and  radiation  heat  transfer,  though  well  understood  and  consistently  modeled  with  high  levels  of 
accuracy,  are  complicated  by  evaporation  (sweating,  respiration)  and  other  phase  change 
processes,  as  well  as  physiological  responses  such  as  vasomotor  reflexes.  Incomplete  anatomical 
data  sets  have  also  prevented  thermoregulatory  models  from  providing  realistic  output. 
Furthermore,  accurate  model  validations  are  difficult  to  perform  because  experiments  are 
invasive  and  boundary  conditions  can  be  difficult  to  quantify.  Animal  models  are  useful,  but  size 
and  morphological  differences  make  it  extremely  difficult  to  extrapolate  from  animals  to 
humans.  Extrapolation  of  experimental  results  from  animal  models  to  humans  is  further 
complicated  by  species  differences  in  thermal  dissipation  mechanisms,  notably  evaporative 
cooling.  The  thermoregulatory  system  of  non-human  primate  species  is  somewhat  similar  to  that 
of  humans,  providing  better  opportunities  for  experimental  studies.  Primate  cooling  physiology, 
particularly  their  ability  to  sweat,  differs  from  that  of  humans  however  [30].  For  this  reason, 
primate  studies,  though  very  useful,  have  somewhat  limited  application  [12,29]. 


Despite  these  factors,  which  thus  far  have  prevented  the  development  of  accurate  anatomical 
human  models,  mathematical  expressions  have  been  derived  which,  though  conceptually  simple, 
accurately  describe  local  temperature  responses  of  homogeneous  tissues  to  thermal  loads. 
Temperature  predictions  are  commonly  based  on  Pennes’  bio-heat  transfer  equation  [11],  which 
is  the  conduction  heat  equation  with  added  terms  to  represent  the  metabolic  heat  load  and 
convective  effect  of  local  blood  flow: 


i,v!r  +  4c,  (r„  -T)+4i=  p,c,  —  (l) 

The  left-hand  side  terms  of  equation  1  represent,  respectively,  conduction  through  tissue  with 
thermal  conductivity  kh  heat  convection  associated  with  blood  flow  rfy ,  and  metabolic  heating. 
The  right-hand  side  term  describes  the  rate  of  energy  increase  of  the  tissue  volume. 

Although  the  theoretical  basis  of  the  Pennes’  model  is  controversial,  it  has  been  validated  in 
numerous  species  and  tissues  and  is  widely  accepted  as  an  appropriate  bio-heat  transfer  equation 
for  most  applications  [2,10].  It  has  also  been  shown  that  blood  flow,  vasodilatation  and 
constriction,  and  other  circulatory  processes  can  be  successfully  modeled  [6],  as  well  as 
metabolic  effects  and  the  role  of  blood  chemistry  [10].  Modem  predictive  techniques  and 
powerful  high-speed  computers  make  application  of  the  bio-heat  equation  feasible  on  animal  and 
human  anatomical  subjects.  Data  sets,  particularly  the  voxelized  data  available  through  the 
Visible  Human  Project  [National  Library  of  Medicine]  and  others  provide  an  anatomical 
description  detailed  enough  to  apply  the  bio-heat  equation  on  a  scale  necessary  to  take  into 
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account  the  effects  of  differing  tissue  properties,  vasomotor  and  sudomotor  responses,  and  the 
convective  effects  of  the  circulatory  system. 


2.0  PHASE  I  TECHNICAL  OBJECTIVES 

The  main  objective  of  the  Phase  I  work  was  to  develop  a  prototype  for  an  anatomically  realistic 
thermal  modeling  tool.  The  prototype,  written  in  C++,  runs  on  both  Unix  and  Windows  systems 
and  is  able  to: 

1 .  Read  voxel-based  digital  anatomical  data  sets  (.raw  files  and  the  corresponding  tissue.txt 
file)  and  produce  a  thermal  nodal  network,  based  on  Pennes’  bio-heat  transfer  equation, 
in  which  each  thermal  node  corresponds  to  one  voxel, 

2.  Use  known  thermal  properties  of  different  tissues  and  materials  to  determine  the  thermal 
capacitance  of  each  node  and  the  conduction  between  nodes, 

3.  Read  and  apply  the  specific  absorption  rate  (SAR),  as  predicted  by  the  FDTD  program,  to 
each  thermal  node  for  a  period  of  time  set  by  the  user, 

4.  Solve  for  the  time  varying  temperature  of  each  node  using  finite  difference  techniques, 

5.  Produce  three  dimensional  temperature  files,  of  the  same  format  as  the  SAR  file,  at  a 
specified  simulation  time,  and, 

6.  Has  been  validated  by  comparison  to  analytical  test  cases  and  experimental 
measurements. 

The  prototype  code  includes  terms  for  metabolic  heating  and  blood  flow  in  the  nodal  energy 
budget,  although  these  will  not  be  available  to  the  user  until  Phase  II.  Our  Phase  I  activities  also 
included  generating  a  database  of  the  thermal  properties  of  different  tissues  and  relevant 
materials.  Thermal  predictions  were  validated  using  analytic  test  cases  and  measured  data, 
obtained  from  the  Directed  Energy  BioEffects  Laboratory  (DEBL)  at  Brooks  AFB,  taken  on  a 
Rhesus  monkey  carcass  subjected  to  far-field  irradiation  under  known  conditions  of  power 
density,  frequency,  and  orientation. 


3.0  TISSUE  PROPERTIES 

Thermal  properties  (thermal  conductivity,  specific  heat,  mass  density,  thermal  diffusivity)  have 
been  identified  for  virtually  all  the  relevant  tissues.  In  some  cases,  there  are  multiple  entries 
corresponding  to  different  species  (rat,  human,  monkey).  [Appendix  A] 

The  metabolic  heat  production  and  blood  flow  per  unit  tissue  volume  are  determined  from  the 
literature,  where  available,  for  most  metabolically  active  tissues.  In  some  cases,  the  metabolic 
heat  production  (per  unit  tissue  volume)  is  estimated  from  O2  consumption  rates. 

For  many  tissues,  blood  flow  rates  are  tabulated  as  a  range  of  values.  This  is  particularly  relevant 
to  superficial  tissues  that  exhibit  significant  vasoconstriction  (or  vasodilatation)  according  to 
thermal  state.  Modeling  this  thermoregulatory  response  will  be  an  important  aspect  of  Phase  II 
work. 


Page  6 


ThermoAnalytics 
Topic:  OSDOO-HPOI 


3-D  Voxel-Based  Bio-Heat  Transfer  Code  SBIR  Phase  I  Final  Report 

Contract:  N00014-01-M-0071 


4.0  PHASE  I  CODE  DEVELOPMENT 

In  this  section  we  will  describe  the  development  of  the  Bio-Heat  Transfer  software  delivered  at 
the  end  of  Phase  I.  Originally,  this  prototype  was  to  be  based  on  ThermoAnalytics’  existing 
thermal  modeling  software,  RadTherm.  Based  on  feedback  from  the  kick-off  meeting,  an 
emphasis  was  placed  on  producing  a  prototype  thermal  solver  that  included  parallelized  code  and 
for  which  the  source  code  could  be  made  widely  available.  For  this  reason,  the  voxel  thermal 
solver  was  implemented  as  a  sub-solution  type  of  our  main  solver  (similar  to  the  tri-diagonal 
solver  we  use  for  ID  terrain  models)  in  such  a  way  that  it  can  also  be  made  standalone. 


4. 1  Model  Import  and  Setup 

The  modeling  process  begins  by  importing  voxel  data  sets  containing  tissue  types  and  SAR 
values.  The  large  size  of  voxel  data  sets  (e.g.,  over  45  million  voxels  for  the  “2mm  man,”  a  293  x 
170  x  939  voxel  representation  of  a  human  body)  requires  that  we  utilize  compact  data 
structures,  avoid  duplicate  data  representations,  and  eliminate  storage  by  recalculating  simple 
quantities  as  they  are  needed.  Fortunately,  the  simple  structure  of  the  voxel  description  lends 
itself  to  this  type  of  optimization.  The  C++  class  (See  Appendix  B.l)  that  holds  voxel  data 
consists  of  variables  that  hold  the  x,  y,  z  dimensions  of  each  voxel,  the  number  of  voxels  in  the  x, 
y,  z  directions  (nx,  %,  nz),  a  one-dimensional  array  dynamically  sized  to  hold  nx  *  %  ♦  nz  instances 
of  a  voxel  data  class,  and  functions  for  accessing  the  voxel  information  at  a  given  index  or 
position.  The  voxel  data  class  contains  storage  for  temperature,  temperature  from  the  last  time 
step,  the  SAR  heating  value,  and  the  tissue  type  of  the  voxel.  This  storage  scheme  is  several 
times  smaller  than  if  voxel  data  were  stored  using  the  data  structures  typically  defined  for  solid 
elements.  (The  storage  for  a  single  solid  element  consists  of  6  integers  that  hold  vertex  indices. 
There  are  roughly  as  many  vertices  as  there  are  solid  elements  and  each  requires  3  floats  to  hold 
the  x,  y,  z  coordinates.) 

Several  parameters  that  control  the  thermal  solution  are  required  and  these  are  read  from  a  text 
file  prior  to  invoking  the  thermal  solution.  These  parameters  are  listed  in  Table  1,  next  page.  A 
sample  text  input  file  is  included  in  Appendix  B.3. 
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Table  1.  Model  Parameter  Input  File 


Parameter 

Units 

Number  of  voxels  in  x,  y,  and  z  dimensions 

n/a 

Voxel  size  in  x,  y,  and  z  dimensions 

mm 

Duration  of  exposure  to  RF  heating 

Seconds 

Duration  of  thermal  simulation 

Seconds 

Power  density  of  RF  heating 

Initial  temperature  of  the  object 

°C 

Time  step  of  thermal  simulation 

Seconds 

Ambient  temperature 

°C 

Surface  convection  coefficient 

W/m2  -K 

4.2  Thermal  Solution 

The  thermal  solver  for  the  Phase  I  prototype  solves  a  variation  of  Pennes’  bio-heat  equation 
which  includes  the  effects  of  local  blood  flow,  metabolic  heating,  convective  heat  transfer,  and 
the  heating  effects  of  RF  energy.  The  following  equation,  repeated  from  Section  1.1  with  added 
terms  <fisar  for  RF  heating  and  <ficonv  for  convective  heat  transfer,  is: 

kyr + fSic .  (r.  -  r)+ 4  +  ft,,  +  ft,  -  P,cn  ^  (2) 

The  spatial  derivatives  in  Equation  2  are  discretized  for  the  voxelized  domain  using  second  order 
finite  difference  expressions,  and  the  Crank-Nicholson  scheme  is  applied  for  discretization  in  the 
time  domain.  The  discretization  yields  the  following  equation: 
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,k 

^ Ph 

2 

~r 

2 

i 

h.Ar 


AxAyAz 

r1,.,  ■ 


7t  «+1  rpn+]  nn  w  n~>n 

oo  —l  ij\k  1  oo  —  1 


n+\ 


-  p£p, 


r«+l  T>n 

ij,k  ~~  1  ijjk 

At 


(3) 


The  superscripts  n  and  n+1  indicate  time  step  n  and  time  step  n+1,  respectively.  The  parameter 
As  is  the  surface  area  of  a  voxel  face  that  is  exposed  to  ambient  air,  and  hs  is  the  coefficient  for 

convective  heat  transfer  between  exposed  voxel  faces  and  the  external  air.  The  parameters  Ax, 
Ay,  and  A z  denote  the  dimensions  of  an  individual  voxel  in  the  x,  y  and  z  directions,  and  At  is  the 
time  step  interval. 

Equation  3  will  be  rearranged  to  solve  for  Tn+'ij,k ,  and  similar  terms  will  be  combined  for  more 
efficient  numerical  evaluation.  The  first  step  is  to  separate  the  terms  that  contain  temperatures 
from  the  n+1  time  step: 
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^1  1^11+]  rpn  + 1  ,  T'/I  +  l  rpn+\  1  I 

— -[1  i-\j,k  - 1  ij,k+l  <+1  JJc-1  ij,k\+ 


-[t'-xj*  -r,,j,k  +TnMJ,k  -T”u*]+ 


Ax 


Ax' 


k/  [r+Vu  -  Tn+'ij,k  +  Tn+',J+hk  -  r+\j,k  ]+ 
[r.j-u  - T",j,k  +  Tni,J+\,k  -Tnij,k\+ 
[r+'u,k- 1  -  Tn+]u,k+Tn+\jMi  -r+'u,k}+ 


Ay 

kt 


Ay 
k , 


Az 

k,  [rpn 


A z1 


[r  -rKJ,k+rUMi  -r, ,,.*]+ 


AcP,  [t;",.,*  -  r"i,Jt  ]+  Acr,  [+  t;,m  -  r  u,  ] 

a  \rr  ,,+'  rpn+\  1  ,  ^sAg  n  rpn  "| 

—  [7W  -1  ij,k J+  A_;;— I/.  -J  u,k J 


h,A 


AxAyAz 
»+i 

"h.jj, 


AxAyAz 

rpn+ 1  nrn 

1  ij\k  1  ij,k 


P, 


At 


Collecting  the  terms  in  Equation  4  that  contain  Tn+ij,k  and  7’",,;,*  yields: 


rw+1 

ijjt 


2kt  2 k,  2k,  «  r~, 

+  ^T  +  ^T  +  r#iLn*  + 


h„  Af  2.p,Cn 

sijjc  Ai.jJt  .  rip 


Ax  Ax  Ay 


-r, 

J,*  - 

i 

K  | 

\rn+\ 

Ax2 

Ly  ‘ 

K 

/7+1 

Ay2 

L 

K  | 

[r«+1, 

Az2 1 

1/  ' 

M  ,  r, 

Ph  AxAyAz  At 

As  2  p.Cn 

Si.jM  sijjc  Pi 


2^+ 2^+ 2^+  + - 

Ax2  Ax2  Ay2  Ph  AxAyAz  At 

[r+\-,,/.*  +r+W*]+- A-Ir '->•>•*  +r<+u,*]+ 

Ax 

[r+l,v-u  +r+l/,7+u]+ A_[r,,y-u  +r,.,+u]+ 
[r+,^-i  +r+I,  ,,,*+.]+ A_[r  + r, ,,,*+.]+ 


Az 


AxAyAz  AxAyAz 

2<r +2#+l,a,,,  =0 

Equation  5  is  rearranged  to  prepare  for  solving  for  Tn+\j,k 


(4) 


(5) 
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Pk  AxAyAz  At 


AxAyAz  At 
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Ax2 

i1  1 

k,  | 
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[rn+ 

A z2  1 

i1  1 

*cp\t: 

m,  r 

Ax 
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+rVu]+ 
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— !y[7’"+,|,y,*-1  +  T"+] tj,k+l]+  ' 2  +  T  ij,*+l]  + 

\z  A z 


V:"V^h :&■•'] 


AxAjyAz 


AxAyAz 


]+ 

1+ 


(6) 


Equation  6  can  be  simplified  by  defining  a  parameter  C,  which  takes  on  several  different  values 
depending  on  the  position  of  the  adjacent  voxels,  the  presence  of  blood  flow,  and  whether  or  not 
any  of  the  voxel  faces  are  exposed  to  convection  with  ambient  air.  The  expressions  for 
evaluating  the  parameter  are  as  follows1 : 


1  The  parameter  C  must  be  modified  for  the  conduction  terms  when  the  adjacent  voxel  has  a  different 
thermal  conductivity  than  the  current  voxel.  In  this  case  the  parameter  takes  the  form  of 
k  k 

C  =  -7 - — - ,  where  k „  is  the  thermal  conductivity  of  the  current  voxel,  and  ka  is  the  thermal 

(*„+U A*2 

conductivity  of  the  adjacent  voxel. 
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k 

C  =  — '-r  { for  conduction  terms  with  subscript  i  ±  1) 

Ax 

k 

C  =  — ‘-j  { for  conduction  terms  with  subscript  j  ±  1) 

A y 

k 

C  =  — {for  conduction  terms  with  subscript  k  ±  1) 

Az" 

C  =  dfC Pi  {for  blood  flow  terms) 
h 

C  =  — {for  convection  on  surfaces  in  the  y-z  plane) 

Ax 

h 

C  =  — -  {for  convection  on  surfaces  in  the  x-z  plane) 

Ay 
h 

C  =  —  {for  convection  on  surfaces  in  the  x-y  plane) 

A z 


(7) 


Applying  the  expressions  in  (7)  to  Equation  6  yields: 

+Ec-7'-"*,+Zc<r<"+ 

1= 1  /= 1  V p ) 

2<ft] mi  jk  +  , 

where  /  is  an  index  for  the  summation  of  all  parameters  C  that  apply  to  voxel  i,j,k.  The  limit  N  is 
the  total  number  of  C  parameters  that  apply  to  voxel  i,j,k,  and  7)  is  the  temperature  associated 
with  parameter  Q,  which  could  be  the  adjacent  voxel  temperature,  ambient  air  temperature,  or 
blood  temperature.  Solving  Equation  8  for  T  n+I iJtk  and  rearranging  some  terms  leads  to  the 
equation  that  is  encoded  in  the  numerical  solution: 


nf1  +  \ 


ij,k 


+ 


/= 1 


2_p£, 

A  / 


rpn 

=  -l  ij,k 


iv 

Zc, 


2  P,C 


Pi 


/« 1 


At 


-»«+! 


UJc  ~ 


J  V 

+ 


2  p,C 


Pi 


At 


+  Y.C,T,"  ~Tnu,k 

/= 1  M 

Yr  2p£p, 

h  1  At 

+ 

> 

(9) 


Equation  9  is  applied  for  each  voxel  in  the  domain  and  the  resulting  set  of  equations  is  solved 
using  an  iterative  successive  over-relaxation  solver.  The  source  code  for  a  portion  of  the  voxel 
thermal  solver  is  included  in  Appendix  B.2. 

As  mentioned  earlier,  the  voxel  thermal  solver  is  implemented  as  a  sub-solution  to  the  main 
thermal  solver  which  is  in  RadTherm.  This  was  done  by  modifying  the  main  thermal  solver  to 
allow  the  addition  of  an  arbitrary  number  of  voxelized  domains.  During  each  iteration  of  the 
main  thermal  solver,  a  single  iteration  of  the  voxel  thermal  solver  is  performed.  This 
implementation  allows  for  the  reuse  of  much  of  the  main  thermal  solver  architecture,  and  will 
also  make  it  easier  to  add  radiative  heat  transfer  boundary  conditions  in  the  future. 
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Thermoregulation  Solver 
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Using  the  solution  parameters  as  specified  in  Table  1,  the  thermal  solver  initializes  each  voxel 
using  the  initial  object  temperature.  The  solution  proceeds  with  solving  for  the  voxel 
temperatures  (as  described  above)  at  each  time  step  for  the  specified  duration.  The  normalized 
RF  heating  values  are  multiplied  by  the  specified  power  density  and  applied  to  each  voxel  from 
the  beginning  of  the  solution  up  to  the  duration  of  exposure.  After  this  time  period  is  over,  the 
RF  heating  is  disabled.  A  constant  convection  boundary  condition  is  applied  at  each  voxel 
surface  that  is  exposed  to  external  air,  using  the  specified  convection  coefficient  and  ambient  air 
temperature.  After  the  thermal  solution  is  complete,  a  dataset  containing  the  temperature  of  each 
voxel  is  written  to  a  binary  file  in  the  same  format  as  the  SAR  heating  values.  Figure  1  contains  a 
flowchart  of  the  solver. 


4.3  Parallelization 

The  objective  for  parallelizing  the  voxel-based  thermal  solver  was  to  be  able  to  solve  large,  i.e., 
memory  intensive,  problems  on  a  distributed  memory  computer  cluster.  We  utilized  the  current 
standard  library  for  writing  distributed  memory  parallel  programs,  the  Message  Passing  Interface 
(MPI).  After  reviewing  the  literature  regarding  the  parallelization  of  various  iterative  solution 
techniques,  we  chose  to  use  Successive  Over-Relaxation  (SOR)  as  our  solution  method, 
primarily  because  it  requires  the  least  amount  of  storage.  All  of  the  solution  methods  considered, 
including  those  that  generally  provide  faster  convergence,  e.g.,  multigrid  and  FFT,  require 
storage  that  is  linear  with  the  number  of  nodes.  However,  the  amount  of  storage  required  per 
node  for  SOR  is  a  fraction  of  that  required  by  the  other  techniques. 

We  followed  the  typical  approach  of  partitioning  the  thermal  nodes  among  the  parallel 
processors  to  minimize  inter-processor  communication.  Since  in  most  test  subjects,  the  x-y  cross 
section  is  the  smallest,  the  problem  domain  was  partitioned  along  the  z-axis.  Each  partition 
contains  one  or  two  extra  x-y  “slices”  of  thermal  nodes,  which  are  a  copy  of  the  nodes  whose 
temperatures  are  calculated  by  the  adjacent  processors).  Since  the  temperatures  of  these  slices 
are  only  updated  at  the  end  of  each  iteration,  the  data  dependencies  across  processor  boundaries 
are  therefore  changed  from  “Gauss-Seidel  style”  (in  which  the  results  from  the  current  iteration 
are  used  as  they  become  available)  to  “Jacobi  style”  (in  which  the  results  from  the  previous 
iteration  are  used  throughout  the  current  iteration).  The  temperatures  of  boundary  nodes  are 
exchanged  between  processors  via  point-to-point  communication.  Figures  2  and  3  summarize  the 
communication  between  processes. 

In  order  to  improve  scalability,  the  data  needed  to  solve  the  entire  thermal  problem  (i.e.,  tissue 
IDs,  temperatures,  and  SAR  values)  are  never  stored  on  a  single  processor.  Instead,  processor  0 
reads  (writes)  buffer-fulls  of  information  that  is  scattered  to  (gathered  from)  the  other  processors. 
In  the  case  of  the  2mm  man  model,  in  which  the  data  is  too  large  to  be  contained  in  a  single 
process  (under  RedHat  7.1),  the  parallel  version  of  the  code  can  be  used  to  run  the  model  in  two 
processes  running  on  a  single  CPU. 

We  installed  the  LAM  6.5.3  version  (free  Notre  Dame  version)  of  MPI  on  our  network  of  PCs 
running  Linux  (5  Intel  P-III  lOOOMhz  PCs  with  512Mb  to  1Gb  of  physical  memory  each  and  2 
AMD  1333Mhz  PCs  with  784Mb  of  physical  memory).  Typical  performance  for  the  serial  and 
parallel  version  of  the  code  is  summarized  in  Table  2  and  Figure  4. 
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Processor  Rank  =  0 


Called  from  main 


Solve3DMeshMaster 

Constructor 


< 


Solve3DMeshMaster 
functions  called  from  main 


Solve3DMeshMaster 
functions  called  from  solver 


Solve3DMeshMaster 
functions  called  from  main 


Solve3DMeshMaster 
Destructor 


Called  from  main 


Processor  Rank  >  0 


Called  from  main 


Solve3DMeshSlave 

Constructor 


Solve3DMeshSlave 

functions  called  from  receive  message  loop 


Called  from  main 


Figure  2:  Overview  of  Communication  within  the  Parallel  Version  of  the  Code 
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Processor  Rank  =  0 

From  Solver  Iteration  Loop 


Processor  Rank  >  0 


From  Solver  Time  Stepping  Loop 


Send  Iteration  Info 
(tag  = 

CopyNewToOld) 

r 

Copy  New 
Temperature 
To  Old 

Figure  3:  Parallel  Solver  Communication  during  Temperature  Solution 
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Table  2:  Typical  runtimes  for  30  seconds  of  simulation  using  a  5  second  time  step 

on  a  single  CPU. 


Model 

#  Voxels 

#  Non-Air 
Voxels 

Runtime 

(seconds) 

PIII-1000 

Runtime 

(seconds) 

AMD-1333 

2.5mm  Monkey 

2,281,600 

397,181 

73 

45 

3mm  Man 

13,987,344 

3,875,766 

637 

393 

2mm  Man 

46,771,590 

13,078,893 

2355 

2011* 

*Not  enough  RA 

M  to  hold  entire  model. 

Parallel  Code  Speedup 


■♦—2.5mm  Monkey 
«-3mm  Man 
2mm  Man 


Figure  4:  Typical  speedup  for  the  parallel  version  of  the  code.  The  first  five 
processors  were  PIII-IOOOs;  the  last  two  were  AMD-1 333s.  All  processes  fit  in 

available  RAM. 

Note:  the  2mm  Man  model  was  too  large  to  run  in  one  process;  for  this  model  two  processes 
were  running  on  a  single  CPU. 

5.0  VALIDATION 

The  Phase  I  Code  has  been  tested  against  measured  temperatures  under  controlled  conditions  and 
against  analytically  derived  data.  In  this  section  we  will  describe  the  process  and  the  results  of 
the  validation. 

5.1  Analytical  Comparison 

The  analytical  test  case  involved  immersing  a  homogeneous  isothermal  sphere  into  an  infinite 
water  bath  and  then  allowing  it  to  reach  equilibrium  temperature. 
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Table  3:  Analytical  Test  Case 


Sphere  Properties 

Initial 

Conditions 

Boundary 

Conditions 

p=  3000  kg/m3 
k-20  W/m  K 
cD=  1001.0  j/kg  K 

Tj,  sphere  =  335  °C 
Tco=  20  °C 

convection, 
h  =  6000  W/m2  k 
no  radiation 

The  exact  solution  for  the  transient,  one-dimensional  application  of  the  heat  equation  for  a 
sphere,  in  dimensionless  form,  is; 


0*  =  ]Tc„exp(-£,2Fo) 


1 


Is" 


sin  (£/*) 


where; 

Fo  =  ar/r02 


c  4[sin(£,)-£,cos(£,)] 

24-sin(2£,) 

T  -T 
T  -T 

initial  sphere  <» 

a  is  the  thermal  diffusivity  of  the  sphere  material,  ris  the  simulation  time,  r0  is  the  sphere  radius, 
r*  =  r/r0 ,  Fo  is  the  Fourier  Number,  Tx  is  the  constant  water  bath  temperature  and  the  values  of 
E,n  are  the  eigenvectors  to  the  equation; 


1-4  cot  4„=Bi 

where  Bi  is  the  Biot  number  [34,35]. 

An  inherent  trait  of  voxel  geometry  is  that  objects  are  modeled  in  Cartesian  space.  Spheres  and 
all  other  non-planar  surfaces,  modeled  with  voxels,  must  be  approximated  by  cubes,  as  shown  in 
Figure  5. 
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Figure  5:  Spheres  modeled  in  Cartesian  and  Spherical  coordinates 

A  result  of  this  approximation  is  that  error  is  introduced  into  the  simulation  of  a  non-planar 
surface,  since  the  modeled  surface  area  of  the  volume  is  not  exact.  This  error  causes  a  slightly 
non-uniform  transient  temperature  distribution  to  be  calculated  over  the  surface  of  the  sphere. 

We  compared  data  points  at  the  center  and  at  the  surface  of  the  sphere.  Figure  6  shows  the 
temperature  comparison  at  the  center  of  the  sphere.  The  edge  temperature  plot  shown  in  Figure  7 
is  averaged  over  the  sphere  surface  in  order  to  account  for  the  voxel  geometry  induced  error. 


Sphere  Center  Temperature  Comparison 


Time  (seconds) 

Figure  6:  Sphere  Temperature  Comparison,  sphere  center 
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Sphere  Edge  Temperature  Comparison 


Time  (seconds) 


Figure  7:  Sphere  Temperature  Comparison,  sphere  surface 

The  error  between  the  Phase  I  Code  prediction  and  the  analytical  solution  is  plotted  in  Figure  8. 
Considering  the  large  initial  temperature,  the  rapid  rate  of  temperature  change,  and  the  voxel 
approximation  to  a  sphere,  the  error  is  considered  to  be  negligible. 


Analytical  Comparison,  Error 


- Center  error 

- Surface  error 


Figure  8:  Analytical  Comparison  Error 


5.2  SAR  to  Temperature  Correlation 

The  locations  of  thermal  “hot  spots”  within  the  voxel  models  correspond  to  the  locations  of  the 
peak  specific  absorption  rate  (SAR),  as  predicted  by  the  FDTD  code.  The  rate  of  tissue 
temperature  increase  can  be  predicted  (approximately)  for  short-term  exposures  by  a  simple 
extrapolation.  For  an  adiabatic  system,  with  negligible  thermal  diffusion,  the  local  rate  of 
temperature  increase  is  given  by  the  SAR  [24]  for  a  given  material.  This  extrapolation  yields  the 
maximum  rate  of  local  tissue  temperature  increase  for  a  “hot  spot”.  Actual  temperature  increases 
will  be  less,  due  to  thermal  diffusion. 


The  view  in  Figure  9  shows  a  mid-plane  saggital  section  of  a  Rhesus  monkey  head.  Figure  9(a) 
shows  normalized  SAR  values  for  the  rhesus  monkey  head,  irradiated  in  MEHK  orientation  at 
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800  MHz,  as  calculated  by  the  FDTD  code.  The  local  SAR  value  at  voxel  (56,40,28)  is 
378.2W/kg  (see  detail  of  Figure  9  for  the  voxel  location).  The  voxel  size  is  2.5mm  x  2.5mm  x 
2.5mm.  Figure  9(b)  shows  the  Phase  I  Code  temperature  profile,  without  diffusion,  for  the 
monkey  head  exposed  for  30  seconds  at  a  power  density  of  500  mW/cm2.  Figure  9(c)  shows  the 
temperature  distribution  with  diffusion.  The  initial  temperature  of  the  head  is  32.5  C  for  each 
test. 


O  :  mid-plane  saggital  view,  detail  of  Voxel  (56,40,28) 


Figure  9:  800MHz,  500mW/cm2,  a)  SAR,  b)  temp-no  diffusion,  c)  temp-diffusion 


The  Phase  I  Code  prediction  for  voxel  (56,40,28),  without  diffusion,  is  shown  in  Figure  10. 
Voxel  (56,40,28)  is  a  spinal  region  of  the  monkey  carcass.  This  voxel  reaches  an  equilibrium 
temperature  of  35.3365  C  following  30  seconds  of  SAR  exposure. 


Spine  Hot  Spot,  Voxel  56,40,28 
No  Diffusion  Test  Case 
SAR  =  378.2  W/kg 


- Voxel  56,40,28 

Temperature  Response 


Time  (seconds! 

Figure  10:  Adiabatic,  No-Diffusion  Temperature  Response 
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The  heat  equation  for  the  no-diffusion  test  case  simplifies 

n  dT 
*  p  dt 


where  Q  is  the  local  SAR  value,  m  is  the  voxel  mass  and  cp  is  the  specific  heat  of  the  voxel 
material.  The  solution  to  this  equation  for  voxel  (56,40,28)  gives  a  temperature  increase  of 
2.8365  C,  -  exactly  equal  to  the  Phase  I  Code  prediction.  The  effects  of  diffusion  will  be 
discussed  in  Section  5.3. 


5.3  Qualitative  Comparison 

A  series  of  qualitative  comparisons  using  voxel  geometry  of  a  rat  were  made.  In  each  test  case, 
the  rat  geometry  consists  of  homogeneous  tissue  equivalent  material  (TEM)  with  the  properties 
described  in  Table  4. 

Table  4:  Tissue  Equivalent  Material  properties 


p,  TEM  mass  density 

EOOOkgAn^ 

cp,  TEM  specific  heat 

3,700  J/kgK 

k,  TEM  thermal  conductivity 

0.535  W/mK 

The  purpose  of  these  tests  is  to  show  that  the  Phase  I  Code  provides  reasonable  and  logical 
results  under  a  variety  of  test  cases.  Table  5  describes  the  tests  which  were  conducted.  In  this 
section,  we  discuss  the  results  of  these  tests. 


Table  5:  Qualitative  Test  Cases 


Test 

Exposure 

Frequency 

(MHz) 

Power 

Density 

(mW/cm2) 

Exposure 

Time 

(seconds) 

Boundary 

Conditions 

Test  Description 

1 

2060 

100 

30 

adiabatic 

Compare  SAR  to 

temperature 

distribution 

2 

2060 

10/100 

300/30 

adiabatic 

Compare  effects  of 
equal  overall 
heating  at  different 
power  densities 

3 

2060 

10 

300 

convection/ 

adiabatic 

Comparison  of 
convection  vs 
adiabatic  boundary 
condition 

4 

900/2060 

100 

30 

adiabatic 

Comparison  of 
equal  heating  at 
different 
wavelengths 

All  results  were  obtained  using  a  +HKE  orientation  for  the  transmitter.  Figure  1 1  illustrates  this 
orientation. 
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Figure  11:  +HKE  orientation  for  the  transmitter 

Figure  12  shows  the  results  of  Test  1.  The  purpose  of  Test  1  was  to  verify  that  the  thermal 
response  to  SAR  heating  is  appropriate.  We  would  expect  the  temperature  profile  over  the  3-D 
voxel  geometry  to  approximately  match  the  FDTD  SAR  distribution,  as  Figure  12  shows.  Note 
that  temperature  diffusion,  due  to  thermal  conductivity  and  the  thermal  capacitance  of  the  tissue 
material,  tends  to  smooth  the  temperature  distribution. 

4 

0.0  SAR  (W/kg  /  W/cni:)  2.0 

Figure  12:  Test  1,  SAR  heating  (left),  temperature  response  (right) 

Figure  13  shows  the  results  from  Test  2,  a  comparison  with  differing  power  densities  and 
exposure  times,  but  equal  overall  heating.  Figure  14  features  the  same  comparisons,  shown  in  a 
plane  view  with  the  highest  SAR  values.  The  purpose  of  this  test  was  to  investigate  the  thermal 
effects  of  varying  SAR  intensities  and  exposure  times.  The  test  shows  that  a  higher  power 
density  for  a  shorter  exposure  time  results  in  higher  peak  temperatures  than  a  low  power  density 
extended  over  a  longer  period  of  time,  even  though  the  overall  energy  absorbed  in  each  test  is 
equal. 
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Figure  13:  Test  2,  SAR  distribution,  low  power  density,  high  power  density  exposures 


Figure  14:  Test  2,  SAR  distribution,  low  power  density,  high  power  density  exposures 

Figure  15  shows  the  results  of  Test  3.  This  test  explores  how  the  thermal  response  to  SAR 
heating  differs  between  an  adiabatic  surface  condition  and  a  convection  boundary  condition.  The 
convection  boundary  condition  tends  to  decrease  the  SAR  heating  effects  near  the  surface  of  the 
rat.  SAR  exposures  that  result  in  heating  of  centrally  located  sections  of  the  rat  will  not  be 
initially  affected  by  this  boundary  condition  because  of  the  time  lag  associated  with  energy 
diffusion. 
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Figure  15:  Test  3,  SAR  distribution,  adiabatic  surface,  convective  boundary  condition 


Figure  16  and  Figure  17  show  the  results  of  Test  4.  This  test  investigates  the  effects  of  SAR 
heating  at  differing  wavelengths.  Figure  16  depicts  heating  at  cellular  phone  frequency,  900Mhz, 
while  Figure  17  shows  the  effects  of  heating  at  microwave  oven  frequency,  2060Mhz.  The 
results  of  this  test  are  similar  to  those  of  Test  1,  where  the  contours  of  the  temperature  response, 
though  diffused,  match  the  distribution  of  SAR  heating.  These  images  also  highlight  the  effect 
that  exposure  frequency  has  on  the  distribution  and  intensity  of  heating. 


Figure  16:  Test  4,  SAR  heating  at  900Mhz,  temperature  response 
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0.0  SARJW/kg/W/cm1)  2.0  Temperature  (°C)  20.85 


Figure  17:  Test  4,  SAR  heating  at  2060Mhz,  temperature  response 


5.4  Monkey  Head  Test  Case  Comparison 

Probe  measurements  from  the  experiment  described  in  Table  6  were  compared  against  Phase  I 
Code  predictions  for  a  Rhesus  monkey  carcass. 


Table  6:  Irradiated  Monkey  Test  Case 


Exposure 

Frequency 

(MHz) 

Exposure 

Orientation 

Exposure 

Power 

Density 

(ni\V/cm2) 

Monkey 

Number 

Comparison 

Location 

Probe 

Track 

number 

Exposure 

Duration 

(seconds) 

Test 

Date 

800 

ME-HK 

500 

798Z 

Spinal  Hot 
Spot 

1 

30 

April 

18 

2000 

The  FDTD  code  was  used  to  predict  locations  within  the  carcass  with  high  SAR  heating.  The 
temperature  probe  measurement  locations  shown  in  Figure  18  were  determined  using  the  FDTD 
results. 
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SENSOR  LOCATIONS,  MONKEY  446D 


800  MHz  X-RAY  with  Probe  Locations 


Figure  18:  FDTD  SAR  prediction  and  probe  locations 

It  is  impossible  to  get  exact  agreement  between  predicted  and  measured  thermal  “hot  spots,” 
since  the  voxelized  Rhesus  monkey  model  used  in  this  test  to  estimate  SAR  heating  is  not  a 
model  of  the  same  monkey  that  was  used  to  take  probe  measurements.  Anatomical  differences 
between  the  test  subject  and  the  model  can  be  seen  in  Figure  18. 

The  Track  1  Spinal  Hot  Spot  probe  is  located  in  a  region  where  SAR  heating  is  spread  over  a 
relatively  large  number  of  voxels.  Consequently,  this  is  the  region  where  probe  measurements 
are  likely  to  agree  with  Phase  I  Code  thermal  predictions.  Figure  19  shows  a  plot  of  spinal 
temperature  measurements  and  Phase  I  Code  predictions  for  the  test  case  described  in  Table  6. 
Although  we  do  not  expect  exact  agreement,  we  do  expect  to  capture  the  proper  trend,  which  is 
the  case. 


Spinal  Hot  Spot,  Track  1 
Probe  Measurement  and  Code  Prediction 


Time  (seconds) 

Figure  19:  Probe  measurements  and  Phase  I  Code  prediction 
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Several  factors  affect  the  accuracy  of  the  prediction  for  this  test  case,  including; 

•  Anatomical  differences  between  test  subjects, 

•  Relatively  large  voxel  size  (2.5mm)  of  the  model  which  can  cause  “hot  spots”  to  be 
averaged  out,  and, 

•  Non-uniform  initial  temperature  of  the  irradiated  test  subject. 


We  attribute  the  Phase  I  Code  under-prediction  in  this  test  case  to,  primarily,  the  anatomical 
difference  between  test  subject  and  the  model. 


5.5  Validation  Summary 

The  Phase  I  Phase  I  Code  has  been  tested: 

•  Against  an  analytical  test  case  (sphere  test),  with  near  perfect  agreement; 

•  Against  a  zero-diffusion  analytical  comparison,  with  perfect  agreement; 

•  For  conceptual  validity  under  several  scenarios  using  rat  geometry;  and, 

•  Against  measured  temperatures  in  an  irradiated  Rhesus  monkey  carcass. 

Based  on  the  results  of  these  tests,  we  consider  the  Phase  I  Code  to  be  an  accurate  prediction  tool 
and  a  validated  baseline  upon  which  the  more  sophisticated  Phase  II  thermal  process  modeling 
will  be  built. 
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Appendix  A.  Tissue  Properties 


iTissue  Type 


Air  (external) 


Air  (internal) 


Species  |kg/m31  |J/kgK] 


Human 


Bone  marrow  (red) 


Bone  marrow  (yellow) 


Brain  cerebellum 


Brain  (whole) 


Brain  (whole) 


Cartilage 


Eye  (cornea)  _ 


Eye(lens) 


Eye  (retina) _ 


Eye  (sclera) 


Eye  (aqueous  humor) 


Gall  Bladder 


Glands 


Heart 


Intestine  (Large) 


Intestine  (small) 


Kidney 


[Human 


(Human 


[Human 


Human 
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Kidney 

Rat 

0.151 

Kidney  cortex 

Human 

1049 

0.547 

1.47E-07 

Kidney  Medula 

Human 

1044 

0.54 

1.48E-07 

Kidney  Pelvis 

Human 

1.37E-07 

Ligaments 

Human 

Liver 

Human 

1050 

3600 

60 

Liver 

Rat 

0.572 

0.0212 

Lung 

Human 

1050 

0.518 

0.00045 

Lymph 

Human 

1030 

Mucous  Membrane 

Human 

Muscle 

Human 

1041 

3500 

0.000516667 

1.65E-07 

0.000338 

Muscle 

Rat 

0.0047 

Nails 

Human 

1300 

Nerve 

Human 

1038 

Pancreas 

Human 

1045 

0.345 

Skin/  dermis 

Human 

1116 

3500 

0.545 

0.00395 

1.4E-07 

Spleen 

Human 

1054 

3720 

0.545 

1.44E-07 

Spleen 

Rat 

0.06167 

Stomach 

Human 

1050 

0.527 

60 

Stomach 

Rat 

0.0275 

Testicles 

Human 

1044 

Tooth 

Human 

2165 

1170 

0.45 

1.78E-07 
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Appendix  B.1  Class  Definitions  for  the  Solver 


/*! 

*  \file  threedmesh.h 

*  \brief  Definition  of  ThreeDMesh  class 

* 

*  This  file  contains  the  definition  of  the  ThreeDMesh  class,  which 

*  stores  the  3-D  voxelized  mesh  needed  for  the  thermoregulation 

*  thermal  solver 
*/ 

/* - 

Copyright  <C)  2001  ThermoAnalytics,  Inc. 


#ifndef  TAI_ThreeDMesh_H 
#define  TAI  ThreeDMesh_H 


/* - 

standard  includes 


*/ 


/* - 

non-standard  includes 


iinclude  "base/taisys .h" 

#include  "taicpp/vector .h" 

#include  "dstruct/threedmeshsize . h" 


public 

#defines, 

enums,  etc. 

public 

typedef s, 

class  definitions,  etc. 

/*! 

*  \brief  Data  class  containing  data  for  each  cell. 

* 

*  The  ThreeDMesh  class  contains  a  vector  of  ThreeDCell  classes. 

* 

*  \sa  ThreeDMesh 
*/ 

class  ThreeDCell  { 


public : 

ThreeDCell {) ; 

ThreeDCell (const  ThreeDCell  &); 
-ThreeDCell ( ) ; 


ThreeDCell&  operator= (const 

float  T ( )  const  { 

float  Told ( )  const  { 

int  tissuelndex ( )  const  { 
float  Qsar()  const  { 

void  setT (float  T) 
void  setTold (float  T) 
void  setTissuelndex (int  idx) 
void  setQsar (float  Q) 


ThreeDCell  &)  ; 

return  vT;  } 
return  vTold;  } 
return  vTissuelndex;  } 
return  vQsar;  } 

{  vT  =  T;  } 

{  vTold  =  T;  } 

{  vTissuelndex  =  idx;  } 
{  vQsar  =  Q;  } 
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private : 

float  vT;  // 
float  vTold;  // 
int  vTissuelndex;  // 
float  vQsar;  // 


Temperature  of  cell. 
Temperature  of  cell  at 
The  tissue  the  cell  is 
The  energy  absorbed  by 


last  time  step, 
made  up  of. 

the  cell  from  RF  sources 


} ; 


/*! 

*  \brief  Data  class  for  storing  the  3-D  voxelized  mesh  information. 

* 

*  This  data  class  stores  the  size  of  the  mesh  (number  of  cells  in 

*  x,  y,  z),  the  cell  dimension  (x,  y,  z)  and  a  vector  of  ThreeDCell 

*  data  classes,  (one  instance  of  ThreeDCell  for  each  cell  in  the  mesh) 

* 

*  \sa  ThreeDCell 
*/ 

class  ThreeDMesh{ 
public : 

ThreeDMesh ( ) ; 

ThreeDMesh (ThreeDMeshSize  SmeshSize) ; 

-ThreeDMesh ( ) ; 

int  resize (ThreeDMeshSize  &meshSize) ; 

int  nCellsO  const  {  return  vMeshSize .nCells ( ) ;  } 

int  nCellXO  const  {  return  vMeshSize  .nCellX  ()  ;  } 

int  nCellYO  const  {  return  vMeshSize .  nCellY  ()  ;  } 

int  nCellZO  const  {  return  vMeshSize .  nCellZ  ()  ;  } 

float  cellSizeXO  const  {  return  vMeshSize . cellSizeX () ;  } 

float  cellSizeY ( )  const  {  return  vMeshSize . cellSizeY () ;  } 

float  cellSizeZO  const  {  return  vMeshSize . cellSizeZ () ;  } 

const  ThreeDCellS  operator () (int  i,  int  j,  int  k)  const; 

ThreeDCell&  operator ()( int  i,  int  j,  int  k) ; 

const  ThreeDCell&  operator  ()  (int  idx)  const  {  return  vCell[idx];  } 
ThreeDCell&  operator  ()  (int  idx)  {  return  vCell[idx];  } 

bool  validlndex (int  i,  int  j,  int  k)  const; 

#if def  H AS__N AME S PACES 

friend  std: : ostreams  operator« (std : : ostream  &s, 

const  ThreeDMesh  &data) ; 

#else 

friend  ©streams  operator« (ostream  &s,  const  ThreeDMesh  &data) ; 

#endif 

void  printSizeO; 
private : 

ThreeDMeshSize  vMeshSize;  //  The  size  and  dimensions  of  the  mesh, 
vector  <ThreeDCell>  vCell;  //  Holds  a  ThreeDCell  instance  for  each  cell. 

private : 

/*!  Disabled  copy  constructor  */ 

ThreeDMesh (const  ThreeDMesh  &) ; 

/*!  Disabled  assignment  operator  */ 

ThreeDMesh^  operator= (const  ThreeDMesh  &) ; 
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/* - */ 

/*! 

/*! 

*  \brief  Retrieves  the  cell  at  the  specified  i,j,k  location. 

* 

*  \return  Returns  a  const  reference  to  the  cell. 

*  \sa 
*/ 

inline  const  ThreeDCell&  ThreeDMesh: : operator () (int  i,  int  j,  int  k)  const 

{ 

return  vCell [ (k*nCellY ( )  +  j  ) *nCellX ( ) +i] ; 

} 

/* - */ 

/*  i 

*  \brief  Retrieves  the  cell  at  the  specified  i,j,k  location. 

* 

*  \return  Returns  a  non-const  reference  to  the  cell. 

*  \sa 
*/ 

inline  ThreeDCell&  ThreeDMesh: : operator () (int  i,  int  j,  int  k) 

{ 

return  vCell [ (k*nCellY ( )  +  j  ) *nCellX ( ) +i] ; 

} 

/* - */ 

/*  i 

*  \brief  Checks  if  the  specified  i,  j,  k  indices  are  valid  indices 

*  for  the  mesh. 

■k 

*  \param  i  Index  in  X  direction.  (should  be  >=  0  and  <  nCellXO) 

*  \param  j  Index  in  Y  direction.  (should  be  >=  0  and  <  nCellYO) 

*  \param  k  Index  in  Z  direction.  (should  be  >=  0  and  <  nCellZO) 

* 

*  \return  Returns  true  if  the  indices  are  valid 

*  \sa 
*/ 

inline  bool  ThreeDMesh : :validlndex (int  i,  int  j,  int  k)  const 

{ 

if  (  i  >=  0  &&  i  <  nCellX ( )  && 

j  >=  0  &&  j  <  nCellY ( )  && 

k  >=  0  &&  k  <  nCellZO  )  { 

return  true; 

} 

else  { 

return  false; 


} 

#endif  /*  end  #ifndef  TAI_ThreeDMesh_H  */ 

/* - 

End  of  file  -  don’t  put  any  code  beyond  the  above  #endif 
- */ 
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Appendix  B.2  Solver  Source  Code 

/*! 

*  \file  threedsolver . cpp 

*  \brief  Implementation  of  the  Solve3DMesh  class.  This  class  is 

*  used  for  exercising  the  thermal  solution  of  voxelized 

*  domain. 

* 

*  The  Solve3DMesh  class  is  intended  to  be  used  along  with  the 

*  main  thermal  solver.  The  main  thermal  solver  has  functions 

*  for  adding  multiple  instances  of  the  Solve3DMesh  class,  each 

*  of  which  would  represent  a  different  voxelized  domain.  The 

*  Solve3DMesh  class  only  has  the  ability  to  exercise  a  single 

*  iteration  of  its  relaxation  numerical  solution.  All  of  the 

*  looping  (iteration  and  time  stepping)  is  handled  by  the  main 

*  thermal  solver. 

* 

*  The  Solve3DMesh  class  uses  Crank-Nicholson  and  SOR. 

* 

*/ 

/* - 

Copyright  (C)  2001  ThermoAnalytics,  Inc. 


/* - 

associated  include 

_ */ 

#include  "threedsolver .h" 

/* - 

standard  includes 

_ */ 

#include  <math.h> 

/* - 

non-standard  includes 

_ */ 

tinclude  ”base/basic .h" 
finclude  "dstruct/matldata.h" 


/* - 

private  #defines,  enums,  etc. 


*/ 


/* - 

private  typedefs,  class  definitions,  etc. 


*/ 


/* - 

private  (static)  global  variables  (used  sparingly) 


*/ 


/* - 

private  (static)  function  prototypes 


*/ 


/* - */ 

/*  ! 

*  \brief  Constructor  for  Solve3DMesh 
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*  \param  mesh  reference  to  the  voxel  domain 

*  \param  material  pointer  to  material  database 

* 

*  \sa  ~Solve3DMesh  () 

*/ 

Solve3DMesh: :Solve3DMesh (ThreeDMesh  &mesh,  MaterialData  ^material) 

{ 

vMesh  =  &mesh; 


material->retrieveConductivities (vConductivity) ; 

material->retrievePropertyForAHMaterials (Material : : Specif icHeat, 

vSpecificHeat) ; 

material->retrievePropertyForAHMaterials (Material: : Density,  vDensity) ; 


vDx  -  mesh . cellSizeX ( ) ; 
vDy  =  mesh. cellSizeY ( ) ; 
vDz  =  mesh. cellSizeZ ( )  ; 


vDx2  =  vDx*vDx; 
vDy2  =  vDy*vDy; 
vDz2  -  vDz*vDz; 

vNCellX  =  mesh . nCellX ( ) ; 
vNCellY  =  mesh.nCellY ( ) ; 
vNCellZ  =  mesh . nCellZ ( ) ; 

vTimeStep  =  1.0; 
vRealTime  =  0.0; 
vDuration  =  10.0; 

} 

/* - 

/*! 

*  \brief  Destructor  for  Solve3DMesh 

* 

*  \sa  Solve3DMesh (ThreeDMesh  &mesh,  MaterialData  ^material) 
*/ 

Solve3DMesh : : ~Solve3DMesh ( ) 

{ 

//  nothing  to  do  for  now. 

} 

/* - 

public  function  definitions 


*/ 


*/ 


/* - */ 

/*  i 

*  \brief  Performs  a  single  iteration  of  the  thermal  solution. 

* 

*  \param  maxDeltaT  The  maximum  change  in  cell  temperature  (deg  C) 

*  for  this  iteration,  (computed  by  this  function) 

* 

*  \return  Returns  SUCCESS  upon  successful  completion. 

*  Other  return  values  indicate  a  failure  occurred. 

*  \sa 
*/ 

int  Solve3DMesh: : iterateOnce (float  &maxDeltaT) 

{ 

int  status  =  SUCCESS; 

bool  steadyState  =  false;  // —  not  allowing  steady  state  solution  for  now. 
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ThreeDMesh  &curMesh  =  mesh(); 

const  int  nNeighborCells  =  6; 

ThreeDCell  *neiCellPtrList [nNeighborCells] ; 

const  float  neiCellDimension [nNeighborCells]  = 

{  dx(),  dx  ( ) ,  dy  { )  ,  dy(),  dz(),  dz{)  }; 

const  float  neiCellDimensionSquared [nNeighborCells]  = 

{  dx2  ( )  ,  dx2(),  dy2  { )  ,  dy2(),  dz2<),  dz2()  }; 

float  deltaT;  //  Change  in  temperature  for  a  node  for  this  iteration 

maxDeltaT  =  -1;  /*  reset  the  maximum  observed  change  in  temperature 

*  in  any  node  during  this  iteration.  */ 

float  absDeltaT;  //  Absolute  value  of  temperature  change. 

//  Loop  over  entire  domain  and  solve  for  the  temperature  of  each  cell 
for  (int  k  =  0;  k  <  nCellZ{);  k++)  { 

for  (int  j  =  0;  j  <  nCellYO;  j++)  { 

for  (int  i  =  0;  i  <  nCellXO;  i++)  { 

ThreeDCell  &cell  =  curMesh(i,  j,  k) ; 

int  tissueldx  =  cell . tissuelndex ( ) ; 

static  const  int  externalAir  -  0;  // —  temporary  define 
static  const  int  internalAir  =  -1;  // —  temporary  define 

static  float  hVal  =  10.0;  // —  temporary  value 

static  float  externalAirTemp  =  293.15;  //—  temporary  value 

if  (tissueldx  !=  externalAir)  {  /*  we  don’t  compute  anything  for 

*  external  air  cells.  */ 
float  kCurCell  =  conductivity (tissueldx) ; 

float  condSum  -  0.0; 
float  condTSum  =  0.0; 
float  condToldSum  =  0.0; 

//  Create  an  array  that  contains  pointers  to  each  neighboring  cell, 
int  n; 

for  (n  -  0;  n  <  nNeighborCells;  n++)  neiCellPtrList [n]  -  0; 
if  (  (i+1 )  <  vNCellX  )  neiCellPtrList [0]  =  &curMesh (i+1 ,  j,  k) ; 
if  (  (i-1 )  >=  0  )  neiCellPtrList [1]  =  scurMesh ( i-1 ,  j,  k) ; 

if  (  (j+1)  <  vNCellY  )  neiCellPtrList [2]  -  &curMesh(i,  j+1,  k) ; 

if  (  (j-1)  >=  0  )  neiCellPtrList [3]  =  &curMesh(i,  j-1,  k) ; 

if  (  (k+1)  <  vNCellZ  )  neiCellPtrList [4 ]  =  &curMesh(i,  j,  k+1) ; 

if  (  (k-1 )  >=  0  )  neiCellPtrList [5]  =  &curMesh(i,  j,  k-1) ; 

//  Loop  for  each  neighboring  cell  and  add  conduction  heat  transfer 
//  contributions  from  each  neighbor.  If  there  is  no  neighbor 
//  on  a  cell  face,  add  convection  heat  transfer  with  the  ambient, 
for  (n  =  0;  n  <  nNeighborCells;  n++)  { 

ThreeDCell  *neighborCell  =  neiCellPtrList [n] ; 

if  (neighborCell)  { 

int  neiTissue  =  neighborCell->tissueIndex ( ) ; 
if  (neiTissue  =~  externalAir)  { 

//  Apply  convection  if  the  neighbor  is  external  air. 
float  conductor  =  hVal/neiCellDimension [n] ; 
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condSum  +=  conductor; 

condTSum  +=  conductor*externalAirTemp; 

condToldSum  +=  conductor*externalAirTemp;  /*// —  modify  if 

*  external  air 

*  changes 

*  with  time.  */ 


} 

else  if  (neiTissue  =-  internalAir)  { 

// —  do  we  do  anything  special  with  internal  air? 

} 

else  { 


//  Apply  conduction.  If  the  neighboring  cell  has  a  different 
//  tissue,  an  eqivalent  conductivity  must  be  calculated, 
float  kVal  =  0; 
if  (neiTissue  !=  tissueldx)  { 

float  kNei  =  conductivity (neiTissue) ; 
kVal  =  2*kCurCell*kNei/ (kCurCell+kNei) ; 

} 

else  { 

kVal  =  kCurCell; 

} 

float  conductor  =  kVal/neiCellDimensionSquared [n] ; 


condSum  +=  conductor; 

condTSum  +=  conductor*neighborCell->T ( ) ; 
condToldSum  +-  conductor*neighborCell->Told ( ) ; 

} 

} 

} 


//  Add  contributions  due 
float  mDotBlood  -  0.034; 
float  CpBlood  -  0; 
float  Tblood  =  311.0; 
float  TbloodOld  =  311.0; 


to  heat  transfer  with 
// —  temporary  value. 
// —  temporary  value. 
// —  temporary  value. 
// —  temporary  value. 


the  blood  in  this  cell. 


float  conductor  =  mDotBlood*CpBlood; 


condSum  +=  conductor; 

condTSum  +=  conductor*Tblood; 
condToldSum  +=  conductor*Tblood01d; 

//  Sum  Q 1 s  (sar  and  metabolic  heating).  Need  to  multiply  by 
//  the  density  since  the  Qfs  are  stored  in  W/kg. 
float  rho  =  density (tissueldx) ; 


float  q  =  rho* (cell . Qsar ( )  +  metabolicHeating (tissueldx) ) ; 

//  Evaluate  and  store  new  temperature  for  this  cell, 
float  Tlast  -  cell.T(); 
float  T; 

if  (steadyState) { 

T  =  (condTSum  +  q) /condSum; 

} 

else  { 

float  Cp  =  specificHeat (tissueldx) ; 
float  tranFactor  =  2*rho*Cp/timeStep ( ) ; 
float  Told  =  cell.ToldO; 


T  -  (condTSum  +  condToldSum  - 

Told* (condSum-tranFactor)  +  2.0f*q)/ 
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(condSum+tranFactor ) ; 

} 


deltaT  =  T  -  Tlast; 


T  =  Tlast  +  relax ()*deltaT; 

absDeltaT  =  fabs (deltaT) ; 

maxDeltaT  =  Max (absDeltaT,  maxDeltaT); 


cell . setT (T) ; 

} 

} 

} 

} 


return  status; 

> 


/* _ */ 

/*! 

*  \brief  Writes  the  thermal  results  to  stdout  for  testing/debugging 

*  purposes. 

*/ 

void  Solve3DMesh: : reportResults ( )  const 

{ 

const  ThreeDMesh  ScurMesh  =  mesh(); 

cout  «  "#Real  Time  =  "  «  realTimeO  «  "  (seconds) \n"; 
cout  <<  curMesh; 


/* - */ 

/*! 

*  \brief  Copies  new  temperature  to  old  temperatures 

* 

*  This  function  is  needed  to  move  new  cell  temperatures  to  old  cell 

*  temperatures  after  a  time  step  is  complete.  This  function  is  called 

*  by  the  main  thermal  solver  at  the  appropriate  time. 

* 

*  \sa  iterateOnce ( ) 

*/ 

void  Solve3DMesh: : copyNewToOld ( ) 

{ 

ThreeDMesh  &curMesh  =  mesh(); 
int  idx; 

for  (idx  -  0;  idx  <  curMesh . nCells {) ;  idx++)  { 

ThreeDCell  &cell  =  curMesh ( idx) ; 
cell . setTold (cell .T ( ) ) ; 

} 

} 


/* - 

public  SLOT  function  definitions 


*/ 


/* - 

protected  function  definitions 


*/ 
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/* - 

private  SLOT  function  definitions 


*/ 


/* - 

private  function  definitions 


*/ 


/* - 

End  of  file 


*/ 
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Appendix  B.3  Sample  Input  File 

#  Input  file  for  thermoreg  thermal  solver. 

######################################################## 

#  Voxel  description 

######################################################## 

#  Specify  the  number  of  voxels  in  X,  Y,  and  Z  directions 
NumberOfVoxelsInX :  51 

NumberOfVoxelsInY:  22 
NumberOfVoxelsInZ :  114 

#  Specify  voxel  sizes  in  mm 
VoxelSizelnX:  1.95 
VoxelSizelnY :  1.95 
VoxelSizelnZ :  2.15 


######################################################## 

#  Numerical  solution  parameters 
######################################################## 

#  Specify  the  duration  of  exposure  to  radiation  in  seconds. 

#  (The  radiation  exposure  will  be  applied  from  time  zero  and 

#  will  drop  to  zero  after  the  specified  time  has  elapsed.) 
DurationOfExposure :  30 

#  Specify  the  duration  of  the  analysis  in  seconds. 

Duration:  30 

#  Specify  the  power  density  of  the  radiation  in  milliwatts/cmA2 
PowerDensity :  100 

#  Specify  the  time  step  in  seconds. 

TimeStep:  5 

#  Specify  the  initial  temperature  of  the  object  in  deg  C 
InitialObjectTemperature :  20 

#  Specify  the  ambient  temperature  in  deg  C 
AmbientTemperature:  20 

#  Specify  the  surface  convection  coefficient  in  W/mA2-K 
Surf aceConvectionCoef f icient :  10 

#  Specify  the  solution  relaxation  value,  (must  be  >  0  and  <  2.0) 

#  Typically  a  value  in  the  range  of  1.0  to  1.8  would  be  used. 
SolutionRelaxationValue :  1.5 
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